#include "yocto/utest/run.hpp"

#include "core/cell-dyn.hpp"
#include "core/full-dyn.hpp"
#include "core/ljc.hpp"
#include "core/fcc.hpp"

using namespace FastMD;


YOCTO_UNIT_TEST_IMPL(dyns)
{
	
	const Real epsilon = 1.0;
	const Real sigma   = 1.0;
	LJCutOff   potential(epsilon,sigma,2.5*sigma);
	const bool isothermal = false;
	for( size_t m=1; m <= 8; ++m )
	{
		const U3D    cr(m,m,m);
		const size_t np = 4 * cr.items();
		const Real   length  = __pow( np, __ONE/__THREE );
		Frame        frame_full(3,np);
		Frame        frame_cell(3,np);
		Frame        frame_clMT(3,np);
		
		frame_full.set_box( V3D(length,length,length) );
		frame_cell.set_box( V3D(length,length,length) );
		frame_clMT.set_box( V3D(length,length,length) );
		
		//-- initialize
		for( size_t i=0; i < np; ++i ) 
		{
			const Real mass = frame_full[ Frame::VELOCITIES ][i].data = __ONE + 4*Rand();
			frame_full[ Frame::POSITIONS ][i].data = __ONE / mass;
			
		}
		
		FCC::Build( frame_full[ Frame::POSITIONS ], cr, frame_full.box, frame_full.ibox, 0.1 );
		
		memcpy( &frame_cell[ Frame::POSITIONS  ][0],  &frame_full[ Frame::POSITIONS  ][0], np * sizeof(V3D) );
		memcpy( &frame_cell[ Frame::VELOCITIES ][0],  &frame_full[ Frame::VELOCITIES ][0], np * sizeof(V3D) );
		memcpy( &frame_clMT[ Frame::VELOCITIES ][0],  &frame_full[ Frame::VELOCITIES ][0], np * sizeof(V3D) );
		
		
		FullDynamics  full_dyn;
		CellDynamics  cell_dyn( frame_cell, 2.8 * sigma );
		CellDynamics  clMT_dyn( frame_clMT, 2.8 * sigma );
		
		std::cerr << "#particles=" << np << ", #cells=" << cell_dyn.GetCells().size << std::endl;
		std::cerr << "Boxes: " << frame_full.box << ", " << frame_cell.box << std::endl;
		full_dyn.eval(frame_full, potential, isothermal);
		std::cerr << "Dyn/Full : " << frame_full.V/np << std::endl;
		cell_dyn.eval(frame_cell, potential, isothermal);
		std::cerr << "Dyn/Cell : " << frame_cell.V/np << std::endl;
		clMT_dyn.eval(frame_clMT, potential, isothermal);
		std::cerr << "Dyn/ClMT : " << frame_clMT.V/np << std::endl;
		
		
		const Vector &afull = frame_full[ Frame::ACCELERATIONS ];
		const Vector &acell = frame_cell[ Frame::ACCELERATIONS ];
		const Vector &aclMT = frame_clMT[ Frame::ACCELERATIONS ];
		Real adiff   = 0;
		Real adiffMT = 0;
		for( size_t i=0; i < np; ++i )
		{
			{
				const V3D da( afull[i], acell[i] );
				adiff += da.normsq();
			}
			
			{
				const V3D da( afull[i], aclMT[i] );
				adiffMT += da.normsq();
			}
		}
		adiff   = __sqrt( adiff )/np;
		adiffMT = __sqrt( adiffMT )/np;
		
		std::cerr << "DiffAccel  =" << adiff   << std::endl;
		std::cerr << "DiffAccelMT=" << adiffMT << std::endl;

	}
	
	std::cerr << "sizeof(Real)=" << sizeof(Real) << std::endl;
	
}
YOCTO_UNIT_TEST_DONE()
